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SUMMARY 

The QR algorithm used to solve over-determined systems of linear equations 
was adapted to execute efficiently on the Control Data STAR-100 computer. Us- 
ing the new vectorized algorithm, the STAR-100 computer solved a system of 
250 equations in 50 unknowns in less than 8.5% of the time it took using the 
original scalar version. This paper describes how the scalar program was adapted 
for the STAR-100 and indicates why these adaptations yielded an efficient STAR 
program. Program listings of the old scalar version and of the new vectorized 
SL/1 version are presented in the appendices. Execution times for the two ver- 
sions, applied to the same system of linear equations, are compared. 

INTRODUCTION 

Programs written in standard FORTRAN language for serial computers do 
run on the Control Data STAR-100 computer, but very inefficiently. To take 
advanatage of the architecture and vector-processing capabilities of the 
STAR-100 computer it is necessary to vectorize the algorithms in these programs. 
Frequently one must rearrange the data and computations. This paper describes 
how the QR algorithm to solve over-determined systems of linear equations was 
vectorized and what factors were considered in developing an efficient STAR 
program. 

The vectorized program utilizes SL/1, a high level language developed by 
NASA's Langley Research Center for the STAR-100 computer. SL/1 incorporates 
many features designed to see that programs it compiles take full advantage of 
the star's architecture and capabilities, including half-word storage and 
arithmetic. SL/1 is compatible with FORTRAN in the sense that programs written 
in either language can call subroutines written in the other. In utilizing the 
program presented in this paper, familarity with some of the notations used in 
the SL/1 language will be helpful. 
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General suggestions concerning the adaptations of algorithms for efficient 
use on the STAR may be found in the paper, "Star Adaptation for Two Algorithms 
used on Serial Computer," by Lona M. Howser and Jules J. Lambiotte (see ref. 1). 

ADAPTATION OF QR ALGORITHM TO SOLVE OVER-DETERMINED 
SYSTEMS OF LINEAR EQUATIONS 

The NASA computer mathematics library presently has a subroutine called 

QRASOS, written in FORTRAN for serial computers, to solve an over-determined 

system of linear equations. This subroutine decomposes the matrix A of the 

system AX=B using Householder transformations. (For details of this algorithm 

see ref. 2). To compute these transformations, it uses subroutines SAXPY, SSCAL, 

SCOPY and function SDOT from the Basic Linear Algebra Subroutines (BLAS) . For 

2 

an mxn (m>_n) matrix A it makes n calls to these subroutines and functions to 
solve the given system. Subroutine calls are very expensive on the STAR-100 
computer. 

In SL/1, a matrix can be stored either column-wise or row-wise. Column 
storage means that elements in one column of the matrix are stored as one 
vector (contiguous locations). Similarly, row-storage means that elements in 
one row are stored as one vector (contiguous locations). In vectorizing this 
algorithm both row and column storage of matrices A and B were tried. 

With row-wise storage, reordering of the scalar-version computations is 
required but use of the inner product macro to decompose matrix A is avoided. 

With column-wise storage the computational steps are the same as in the scalar 
versions with vector instructions replacing scalar instructions, but use of the 
dot product macro is required. It was expected that, because of the avoidance 
of the dot product macro, the row-wise approach would offer a considerable 
saving in CPU time. 

Test results show that in using the STAR computer for this algorithm both 
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vectorized versions offer considerable CPU time savings over the scalar program, 
but that contrary to expectations column— wiss storage is more efficient than row- 
wise storage (see table). 


Size of 
Matrix A 

CPU TIME IN SECONDS 
TO DECOMPOSE 

CPU TIME IN SECONDS 
TO DECOMPOSE AND SOLVE 

New Vectorized 
Version (Colunm 
Storage) 

New Vectorized 
Version (Row 
Storage) 

Old Scalar 
Version 

New Vectorized 
Version (Column 
Storage) 

250 X 

200 

2.169 

2.695 

26.34 

2.239 

120 X 

100 

.405 

.588 

3.396 

.433 

250 X 

50 

.158 

.863 

2.223 

.178 

100 X 

10 

.006 

.069 

0.056 

.009 

200 X 

30 

.053 

.416 

0.698 

.063 


Algorithms for both row-wise storage and colirnm storage to decompose the 
matrix A are given. Back substitutions are not discussed here. In both 
algorithms, A is an mxn matrix and WK is a vector of length n. 

COLUMN STORAGE 

When matrix A is stored column-wise, the decomposition of A is achieved 

t'Vi 

as follows: (Note: All references to i column refer to column entries on 

and below the diagonal) . 

(1) Take the inner product of i^^ column with itself and store in WKj^ 

(2) Take the square root of WK^ 

(3) If WK^ = 0, go to step (10) 

(4) WK^ = WK^ X Sign of A^^^ 

(5) Divide column i by WK^ 

(6) Add 1 to A^ ^ 
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These 6 steps compute the Householder transformation for column i. 

To apply this transformation to the columns K = i+l,....n do the follow- 
ing steps; 

(7) Take the inner product of column i with column K and store the 
result in t 

(8) Divide t by A^.^. and then store the negative of the result in t 

(9) Multiply column i by t and add the result to column K 

(10) Store A. . in t, -WK. in A. and t in WK^ 

When i=n perform steps 1 thru 6 and 10. 

ROW STORAGE 

When the matrix A is stored row-wise, the decomposition of A is achieved 

as follows: (Note: In the steps 1 thru 6 below, all references to the row j 

til 

in the i step of decomposition refer to the entries on and to the right 

of the diagonal. All references to the vector WK refer to its i^^, (i+1)^^ 

....n elements. In steps 7 thru 10 all references to the row j in the 

i^^ step of decomposition refer to the entries to the right of the diagonal 

th til til 

and all references to the vector WK refer to its (i+1) , (i+2) ••••n elements), 

th 

At i^ step of decomposition (i=l,2, . . . .n-1) . 

(1) Set WK=0 

(2) For j=l,2....m, multiply row j by A-i . and add the result to WK 

(3) Take the square root of WK^ 

(4) If WK^=0 go to step 11 

(5) Multiply WK. by sign of A. . 

1 1 » 1 

(6) Divide A, . by WK. and add 1 to the result 

1,1 1 

(7) Divide WK by WK^ and add row i of A to WK 

(8) Divide WK by -A. . 

(9) For j=l+l, i+2,....m 

Divide A . . by WK. 

1,1 1 

■■ . ..... . 



(10) For j=i, multiply WK by A, ^ and add the result to row 

j of A 

(11) Store A^ ^ in t, -WK in A^^^ and t in WK^ 

When i=n, perform Steps 1 thru 6 and 11. 

WHY ROW-STORAGE IS SLOWER THAN COLUMN- STORAGE 
As pointed out earlier, if the matrix A is stored row-wise, the use of the 
inner product macro is avoided and the computation of the Householder trans- 
formations and their application to other columns at each step of the decomposi- 
tion is accomplished by the use of a vector multiplication by a scalar and then 
a vector addition. This should result in a considerable savings of the CPU time 
for a large matrx. However, our numerical experiments show just the opposite. 

This can be explained as follows: When an mxn (m>=n) matrix A is stored row- 

wise, the vector lengths in that algorithm are proportional to n, the smaller 
dimension. On the other hand, for column-wise storage the vector lengths are 
proportional to m, the larger dimension. Equivalently, we see that the row- 
stored algorithm requires more vector start-ups ((m-n) (m-n+I)/2 more) to do 
the same number of total computations as the column-stored algorithm, thus requir-.. 
ing more CPU time to do the same amount of work. 

Another factor which makes the row-stored algorithm slower is that the 
transformation elements are stored in the columns of the decomposed matrix. 

If the matrix is stored row-wise, this leads to additional scalar computations, 
notably in step 9 of the algorithm. This slows down the computations considerably. 
Also, if m is large, then not all m vectors in row-wise storage reside in the 
memory at the same time. Because of need to reference different columns at 
different steps of algorithm, this could lead to excessive paging. Thus, any 
advantage gained by avoiding the use of the inner product in the row— wise storage 
is offset by the need to perform many scalar operations, more iterations and 


excessive paging. 
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APPENDIX A 


SL/1 Coding of QR Algorithm 



MODULE Ml » iOPT-l^tSCURCECl > 

/X 

/^PURPOSE 


/* 

/X 

/m 

/X 
/X 
/X 

/XUSE 
/X 
/X 
/X 

/XPARAMETERS 

/X 


TO SOLVE M SIMULTANEOUS EQUATIONS IN H UNKNOWNS WITH IP X/ 

RIGHT HAND SIDES SO THAT ^HE SOLUTIONS ARE THE EEST POSSIBLE^/ 
FIT IN THE LEAST SQUARES SENSE. THE ROUTINE USES HOUSE- 
HOLDER TRAHSr CRMATIONS TO PERFORM. ^hE OR DEC0MP05 IT ! CH 
OF THE COEFFICIENT MATRIX. 


CALL Q4aRASOS(MAXM,MAXN,M,N,IP,A,B,WT, J09,X,R3D,SUM,WK, lERR) 


/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 

/X 


MAXM 


MAXN 


M 


IP 


AN INPUT INTEGER SPECIFYING THE FIRST DIMENSION OF THE 
A,B, AND RSD ARRAYS IN THE CALLING PROGRAM. MAXM MUST 
BE GREATER THAN OR EQUAL TO M. 

AN INPUT INTEGER SPECIFYING THE FIRST DIMENSION OF THE 
X ARRAY IN THE CALLING PROGRAM. MAXN MUST BE GREATER 
THAN OR EQUAL TO N. 

AN INPUT INTEGER SPECIFYING THE HUMBER OF ROWS OF THE 
A AND S ARRAYS. M MUST BE GREATER THAN OR EQUAL TO N. 

AN INPUT INTEGER SPECIFYING THE NUMBER OF COLUMNS OF 
THE A ARRAY. 

AH INPUT INTEGER SPECIFYING THE HUMBER OF COLUMNS OF 
THE B ARRAY. 

AN INPUT/OUTPUT TWO-DIMENSIONAL ARRAY WITH FIRST OIMEM- 
SICH EQUAL TO ilAXM AND SECOND DIMENSION AT LEAST N. 

OH INPUT, A MUST CONTAIN THE MATRIX OF COEFFICIENTS OF 
THE SYSTEM OF EQUATIONS. ON OUTPUT, A CONTAINS INFOR- 
MATION DESCRIBING THE QR DECOMPOSITION OF A. 

AN INPUT TWO-DIMENSIONAL ARRAY WITH FIRST DIMENSION 
EQUAL TO MAXM AND SECOND DIMENSION AT LEAST IP. 

THE COLUMNS OF 0 MUST CONTAIN THE IP RIGHT HAND SIDE 
VECTORS. 


WT 


X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

V/ 

X/ 

X/ 

xy 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 


AH INPUT OUE-DIMENSIONAL ARRAY OF WEIGHTS. IT MUST 
HAVE LENGTH AT LEAST M. IF WEIGHTING IS DESIRED, 

THE FIRST M LOCATIONS MUST CONTAIN REAL NUMBERS GREATER '«/ 
THAN ZERO. IF WEIGHTING IS NOT DESIRED, WT C13 MUST BE X/ 

X/ 
X/ 
X/ 


A NEGATIVE REAL NUMBER. 

JOB AN INPUT INTEGER SPECIFYING RESULTS TO BE COMPUTED. 


- 1 COMPUTE SOLUT ! ONS ONLY . 

-E COMPUTE RESIDUALS ONLY. 

«3 COMPUTE BOTH SOLUTIONS AND RESIDUALS. 

AH OUTPUT TWO-DIMENSIONAL ARRAY CONTAINING THE SOLU- 
TIONS. X MUST BE DII-1EHSIOHED WITH FIRST DIMENSION 
EQUAL TO MAXN AND SECOND DIMEHSICH AT LEAST IP. IF 
SOLUTIONS ARE DESIRED INTO MATRIX B THEN MAXH MUST BE 
EQUAL TO MAXM. FOR THIS PARTICULAR CASE. 


X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 



/* 

R3D 

AH OUTPUT TWG-uIHLW5i OHhL ARF;AY CGWTmIKING THE RE3ID- 

SJ/ 

/* 

* 

UALS. R5D i*iU3T BE DI HENS I OKED WITH FIRST DIHEHSIOH 

Sc/ 

/ya 


EQUAL TO HAXM AND SECOND DIMENSION AT LEAST IP. 

*/ 

/'ik 



*/ 

/a 

3UH 

AH OUTPUT ONE -DIMENSIONAL ARRAY CGNTA IKING THE WEIGHTED 

Si/ 

/* 


dU'*i3 ut ^viUMKC-i* ut t HE RcoIuLALx^. i>uM i*»UST BE DIMEn— 

*/ 

/St 


SIGHED AT LEAST IP. 

St/ 

/» 



Sc/ 

/* 

WK 

A GHE-DIMEKSIuKhL WORK ARRAY- WHICH MUST BE DiHEKSIOHED 

St/ 

/3r 


AT LEAST H. OH OUTPUT, WK COKTAiHS IKFORMATIOH OH THE 

Si/ 

/Sc 


GR DECui'iPuSITIOK OF A. ^ 

Sc/ 

/* 



Sc/ 

/W. 

lERR 

AH INTEGER ERROR CODE. 

Si/ 

/Sc 



Sc/ 

/Sc 


•0 KO ERROR DETECTED. 

Sc/ 

/St 


•1 H IS GREATER THAN M. 

*/ 

/Sc 


•2 THE DECOMPOSED MATRIX IS SINGULAR. 

St/ 

/St- 


■3 WEIGHTING WAS REQUESTED AND ONE OR MORE WEIGHTS 

Sc/ 

/St 


IS KEGATiVE. 

Sc/ 

/* 



Sc/ 

/St 



Sc/ 

/St 

SOURCE HAi*1PTOK INSTITUTE, HAMPTON VA. 

Sc/ 

/* 



St/ 

/* 

LANGUAGE SL/i. 

«/ 

/St 



St/ 

/Sc 

DATE RELEASED JANUARY 18,1360 • 

»/ 

/St 



ifc/ 

/St 



St/ 



EHTRY PROCEDURE unGRASGS CMhaH ,MhXH ,H,H , IP ,A ,B >:mT , JOB ,X ,R3D , 

SuM,wK > I ERRO s 

REAL VECTOR [imAaMj AkPAYvK; A; 

REAL VECTOR iHAXKj ARRAYCIF} Xj 
REAL VECTOR lMAXi*ij ARRAYCIP5 G.R5D) 

REAL VECTOR IMl wT ; 

REAL VECTOR iHl WKj 
REAL VECTOR iIP3 SUM; 

AUTOKiATiC REAL T% 

INTEGER I , J ,K ,L ,M,H , I F ,MAXM, I ERR ,i tAXN , JOB j 


/St 




/St 

CHECK 

F3R M 

LESS 

/Sc 

IF 

M < H 

THEN 


/* 
/* 
/ * 

/* 

/5K 

/*■* 


GO TO LABI 

ELSE 

CKECR FOR WEIGHTING 
IF WTtii >"0 THEN 
CHECK FOR ILLEGAL WEIGHTS 


5fc’/ 

5K/ 


5ic/ 

jfc/ 

3K/ 


I : ” SELLT CwT , OD ; 

IF I < ?i THEN IERR:«3; 

GO TO LABi 

ELSE 

WT r i f HI » •= SORT CWT t i : M 3 j ^ 

FOR TO H DO 

A Cl } t i ; M 3 : *■ A C 1 > L i : i**i j i»:WT [ i ; M J 5 
EHuF % 

FOR i : “ i TO IP DO 

B C iy r i : I'l 3 ; = B C W r j. : IM j iiiWT I I : !*t 3 5 

THFiF • 
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t.nu i ; 
EHDi$ 


EHDI 


/5K CALL &4SGRDC TO DECuMruSE MATRIX A. 

/* 

CALi. Q4 SQRDC Ch ,MAaM ,M , H , wK3 j 

/3K 


3 ^/ 


/5K 

/* 

/» 


/* 
/ * 
/* 


/ ; 
/jfc 


*/ 


/* CALL Q45QR5L TO SOLVE IP RIGHT HAHu SIDES. 

/« 

CALL u*ISQRSLCHMXM,HAXN,M,H,IF,M,B/iaT,JoB,K,R30,SUM,WlC, IERR>, 

IF IERR>0 THEH IERR:«2 
EHDi ) 

LABlt EKDPj 

/M «^«a^*^***/ 

*/ 

PROCEDURE QnSuRDC CA ,MAXM,im,K ,WR} j 

REAL VECTOR CMAXMi ARRAYCH> Aj 
REAL VECTOR CH3 MX; 

AUTOMATIC REAL Tj 

^ HTtGtK I }U,K»L >M,H,IP , MAXM ^ I ERR ,MAXH , JOB 5 


COMPUTE KH TRAHSFuRMATIOrf FOR COLUMK I 

FOR i:*l TO K“i DO 

WK[Ijx» ACI5ti:Mj .DOT. ACIjiIiMjj 
WK i I 1 : «■ SORT CWK t I j 3 5 
IF MX i I j > 0 THEH 

WKtl j 2 « WKE I j^ABSCACw L I j>/ACI3 E I 3 5 
A Cl > E I : M 3 : *■ A C I } E I ; M 3 /MK E I 3 j 
ACIjEI3i*ACIjEI3'»'i. j 

Hr'rLY tiH 1 KHrts#r ukMm i luH lO kEST OF THE CuLUriHS 


*/ 

a/ 


/ 

*/ 




J t “ i ♦ i 5 

FOR Kx* J TO H DO 

T:* HCi3*CI:M3 .DOT. ACK/EIxMI; 
T:-'T/ACT^EI3j 

hCX^EI:Mj THt ACijlJxMjj 

EHDF ? 

EHDI } 

EHDF 5 

WX 1 H 3 2 • A C H ^ E H 2 1'1 3 • DOT . A C H j E H : M j ; 
wXEHjs** SQRTCWXE H j 5 
IF WX E H 3 > 0 THEH 

WX E H 3 2 * WX E H AB5 CACHjEH3}/ACH^EH3j 
A C rO L K : M J 2 A C H> E H 2 M j / WX E H 3 5 
ACH3EHjx“ hCH!)ENj-> 1 j 

EHDI J 

EHDF 5 * 

PROCEDURE Q43QRSL CilAXM , M mXH , M , K , I P , A , B , W7 , JOB , X , R3D , 

SUiM , WX , I ERR jr 5 

REAL VECTOR Ei‘iAXM3 ARRAY CHj A; 

REAL VECTOR Ei'jAaH 3 ARRAY CIF3 Xj 
KC.AL vuCTOR EMAXi‘i3 ARRAY Cir) B,R3Ds 
REAL VECTOR IW3 WX 5 
REAL VECTOR EM3 WT ; 

REAL VECTOR LiP3 SUMj 

i H I tCt^R I , J , X , L , 1 1 > fi , I p ,?mA aM > I ERR ,HAXH , JOB j 
AUTOMATIC REAL Tj 
J ERR 2 -0 J 
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/« 

/3K 

/* 


•SPECIAL ACTION WHEN 
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/M 

/X 

/3K 


/* 

/3K 

/3K 


IF M - 1 THEN 

IF WKIll - 0 THEN 
IERR;-1-, 

GO TO LAB4 
ENDI, 

IF JOB 02 THEN 

FOR I:-l TO IP DO 

XCIM13x“ BCn£lj/ACl)mi 
ENDFj 
ENDU 

IF JOB <> 1 THEN 

RSDCi:>Cl:IP]:-3.0j 

ENDIs 

GO TO LA04-, 

ENDI, 


COMPUTE TRANS CQ^«B 


K/ 

5K/ 

*/ 


FOR I I -1 TO H DO 
IF wvcm <> C THEN 
FOR J;-l TO IP DC; 

Ti- ACmi:M3 .DOT. 

Tt- -T/ACIDin 5 

BCJ)II:M3«- ♦ T^ACnCItMIj 

ENDF; 

ENDI-, 

EHDFj 

FOR I j-1 TO IP 00 

XCnil:N3 :-BCntl»H3 > 

ENDF j 

IF JOB > 1 THEN 

COMPUTE THE RESIDUES 

FOR I v»i TO IP on 

RSDCi:[i:H3 s*O.Oj 
EHDFs 

FOR I s«l TO IP DO 

PSDCI)CK:M1 :«BCnCKjM3 

ENDFt 

FOR K :« H DGWNTO 1 DO 

FOR L :-l TO IP DO 

T»«ACKD CK;M3 .DOT. RSDCL) IK:M3 t 
IF WKIK3-0 THEN 

1ERR:-K*> 00 TO LAO 4 

END I 5 

T:--T/ACKDIK3 5 

RSDCL>IK:M3:-RSDCL>[K;M3 + T'«>A OO £ K ; M3 5 
SUMCL3 * -RSDCL) 1 1 :M3 . DOT. RSDCLj C 1 :M: ; 
EHOFj 

EMDF-, 

IF WTr.3 > 0 THEN 

FOR I 1 TO IP DO 

RSD C m 1 : M3 : - FED C I> [ 1 : M 3 /WT [ 1 : M3 5 
ENDF 5 

WTCl :M3 : «WTC 1 :M3v'WTt 1 :M1 5 
ENDI-, 


3K/ 

*/ 

3K/ 


EUDI i 



IF JOB 02 '^HEH 


/* »/ 

/* COMPUTE THE SOLUTIOHS */ 12 

/'< «/ 

FOR I:- H DOWHTO 2 DO 
IF WK[ I }-0 THEM 

lERRi-I 5 GO TO LAB4 
ELSE 

FOR J»- 1 TC IP DO 

xcj)cn:- -xcj^m/wKm 5 

Ti« -xcvomj 

XCJ)lliKli- XCJ)C1»K] + T3KftCntl:K3} 

ENDF j 
EHDIi 
EHDFs 

FOR I:-l TO I? DO 
IF WKm-0 THEN 

lERRs-l; GO TO LAB4 
ELSE 

xcmiis- -xcniii/wKLi] j 

END I ^ 

ENDF 5 

ENDI j 


/* 




/* 


SAVE THE TRANSFORMATION 

*/ 

/* 

LAB4s 

FOR 1 : -t TO N DO 

*/ 


Tr-ACUtlls ACnm:*“i:KCI35 WKt!3:*T? 
ENDFj 

EHDPj 

ENDMj 

wm w «» 
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APPENDIX B 


FORTRAN Coding of QR Algorithm 



SUBROUTINE QRftSOSCMAXil ,MAXN , M , N , IP , ft ,C ,WT , JOB ,X , R9D , SUM ,WK , lERRD QRASOOlO 


C* 

C3K PURPOSE 
C» 

c* 

C3K 
C5K 
C* 

C* 

C* USE 


AQRAS003D 
*QRAS0040 

TO SOLVE M SIMULTANEOUS EQUATIONS IN N UNKNOWNS WITH IP *QRAS0050 

RIGHT HAND SIDES 50 THAT THE SOLUTIONS ARE THE BEST POSS I BLE»QRA30060 
FIT IN THE LEAST SQUARES SENSE. THE ROUTINE USES HOUSE- *QRAS0070 
HOLDER TRANSFORMATIONS TO PERFORM THE QP DECOMPOSITION 
OF THE COEFFICIENT MATRIX. 


^DRASOOBO 
*QRAS0090 
«QRAS0100 
XQRASOllO 
9KQRAS0120 
*QRAS0130 
»QRAS0140 
»QRAS0150 
^QRASOISO 

MAXM AN INPUT INTEGER SPECIFYING THE FIRST DIMENSION OF THE 5KQRAS0170 
A.B, AND RSD ARRAYS IN THE CALLING PROGRAM. MAXM MUST ^QRASOIBO 
BE GREATER THAN OR EQUAL TO M. ?*^QRAS0190 

3«QRAS0200 

MAXM AN INPUT INTEGER SPECIFYING THE FIRST DIMENSION OF THE 3KORAS0210 
X ARRAY IN THE CALLING PROGRAM. MAXN MUST BE GREATER 3KQRAS0220 
THAN OR EQUAL TO N. *QRAS0230 

V.QRA50240 

M AN INPUT II^TEGER SPECIFYING THE NUMBER OF ROWS OF THE 5»^QRAS0250 
A AMD B ARRAYS. M MUST BE GREATER THAN OR EQUAL TO N. ;KQRAB0260 

^QRAS0270 
5KQRAS0280 
5KQRAS0290 
'i=QRAS0300 


C* 

C* CALL QRASOSCMAXM,MAXN ,M,H , I P , A ,B ,WT , JOB ,X ,RSD ,SUM,WK , lERR) 
C» 

C5T PARAMETERS 
C« 

C3K 

C3K 

C» 

c* 

C3T. 

c* 

CX 
C5»5 

cy. 

C2K 

CX 
CX 

cw. 


N AH INPUT INTEGER SPECIFYING THE NUMBER Or COLUMNS OF 
THE A ARRAY. 

IF AM INPUT INTEGER SPECIFYING THE HUMOER OF COLUMHS OF 




^ORA50310 



Cx 

Cx 

Cx 

CX 

Cx 

C* 

c* 

Cx 

Cx 

CX 

CX 

CX 

CX 

cx 

CX 

CX 

CX 

CX 

cx 

CX 

cx 

CX 

cx 

cx 

cx 

cx 

cx 

CX 

CX 

CX 

Cx 


THE B ARRAY. XQRASC320 

XQRAS0330 

A AH IHPUT/OUTPUT TWO- 0 IMEHS I GHAL APRAY WITH FIRST OIMEH- XQR,AS0340 
SION EQUAL TO MAXM AND SECOND DIMENSIOM AT LEAST N. XQRAS0350 

ON INPUT, A MUST CONTAIN THE MATRIX OF COEFFICIENTS OF XQRAS0360 
THE SYSTEM OF EQUATIONS. CH OUTPUT, A CONTAINS INFOR- XQRAS0370 
MATIOH DESCRIBING THE QP DECOMPOSITION OF A. XQRAS0380 

XQRAS0390 

B AH INPUT TWO-DIMENSIONAL ARRAY WITH FIRST DIMENSION XQRAS0400 

EQUAL TO MAXM AND SECOND DIMENSION AT LEAST IP. XQRAS0410 

THE COLUMNS OF 9 MUST CONTAIN THE IP RIGHT HAND SIDE XQRAS0420 
VECTORS. XQRAS0430 

XQRASG440 

WT AN INPUT OHF.-DIMEMSIOHAL ARRAY OF WEIGHTS. IF WEIGHT- XQRAS0450 

ING IS DESIRED, WT MUST HAVE LENGTH AT LEAST M, AND XGRAS0460 
THE FIRST M LOCATIONS MUST CONTAIN REAL NUMBERS GREATER XQRAS0470 
THAN ZERO. IF WEIGHTING IS HOT DESIRED, WT CAH CONSIST XQRAS0480 
OF A SINGLE LOCATION WHICH MUST CONTAIN A NEGATIVE REAL XQRA3049C 
HUMBER. XQRA50500 

XQRAS0510 

JOB AN INPUT INTEGER SFECIFYIHG RESULTS TO BE COMPUTED. XCRASOSZO 

XQRAS0530 

-1 COMPUTE SOLUTIONS ONLY. XQRA50540 

»2 COMPUTE RESIDUALS ONLY. XQRAS0550 

-3 COMPUTE BOTH SOLUTIONS AND RESIDUALS. XQRAS056C 

XQRA50570 

X AN OUTPUT TWO-DIMENSIONAL ARRAY CONTAINING THE SOLU- XQRAS058G 

TIGHS. IF JOB=l OR J0B«3, X MUST BE DIMENSIONED WITH XORAS0590 
FIRST DIMEHSICM EQUAL TO MAXN AND SECOND DIMENSION xQRASOGCO 

AT LEAST IP. IF JCB«2, X CAN BE A DUMMY PARAMETER. XQRAS061C 


XQRAS0B20 



0^ 

C5K 

c* 

cx 

CX 

c* 

cx 

cx 

cx 

Cx 

CX 

CX 

CX 

CX 


RSO AN OUTPUT TWO-DIMENSIONAL ARRAY CONTAINING THE RESID- 
UALS. IF JOB-2 OR J09-3, RSD MUST BE DIMENSIONED WITH 
FIRST DIMEMSION EQUAL “I’D MAXM AMD SECOND DIMENSION 
AT LEAST IP. IF JOB-1, RSD CAN BE A DUMMY PARAMETER. 

SUM AN OUTPUT ONE-DIMENSIONAL ARRAY CONTAINING THE WEIGHTED 
SUMS OF SQUARES OF THE RESIDUALS. IF J09-2 OR JOB-3, 
SUM MUST BE DIMENSIONED AT LEAST IP. IF JOB-1 . SUM 
CAN BE A DUMMY PARAMETER. 

WK A CHE-DIMENSIONAL WORK ARRAY WHICH MUST BE DIMENSIONED 
AT LEAST N. ON OUTPUT, WK CONTAINS INFORMATION OH THE 
QR DECOMPOSITION OF A. 

lERR AN INTEGER ERROR CODE. 


XQRAS0630 

X0RAS0640 

XQRAS0550 

XQRAS0660 

XQRAS0670 

XQRAS0680 

XQRAS0690 

XQRAS0700 

XQRAS0710 

XQRAS0720 

XQRAS0730 

XQRAS0740 

XQRAS0750 

XQRAS0760 

XQRAS0770 

XQRAS0780 


CX 

CX 

cx 

CX 

CX 

CX 

CX 

cx 

CX 

Cx 

CX 

CX 

CX 



-0 

NO 

ERROR DETECTED. 

XGRAS0790 


-1 

N I 

S GREATER THAN M. 

XQRAS0800 


-2 

THE 

DECOMPOSED MATRIX IS SINGULAR. 

XQRAS0810 


-3 

WEIGHTING WAS REQUESTED AND ONE OP MORE WEIGHTS 

xQRA5O02a 



IS 

NEGATIVE. 

XQRAS0830 





XQRAS0S40 

REQUIRED 

ROUTINES 

NORMS ,SQRDC2 ,SQRSL2,SAXPY1 ,3D0T1 ,SSCAL 

XQRAS0850 




SCOPY 

XQRAS0860 





XQRAS0870 

FORTRAN 

FUNCTIONS 

A8S,AMAX1 ,MIH0 ,MOD ,S IGN ,SQRT 

XQRAS0830 





XGRAS0890 

SOURCE 



COMPUTER SCIENCES CORPORATION, 

XQRAS0900 




HAMPTON, VA. 

XQRA50910 


cx 

; CX LANGUAGE 


FORTRAN 


XQRA50920 

XQRAS0930 


i*:nD6qnQJ A 



DATE RELEASED 


AUOUET 1, 1978 


:kQRAS0950 


cm 


cm 


mORASOSeO 


cm LATEST REVISION OCTOBER 10, 1978 


iKQRAS097 0 


cm 


XQRA50980 

CM 


XQRAS0990 


--aRASlOOO 


DIMEHSlOH ACMAXM,1> ,BCMAXM,15 ,XCMAXN,1D ,RSO CMAXM , ID ,WTC1D ,MKC1D 

QRASlOlO 


DIMENSION SUMCID 

QRAS1020 


lERR - 0 

QRAS1030 

0 


31RAS1040 

c 


QRAS1050 

c 

CHECK FOR M LESS THAN H. 

QRAS1060 

c 


QRAS1070 


IFCM .GE. ND GO TO 10 

QRAS1080 


lERR - 1 

QRAS1090 


GO TO 160 

QRASllOO 

c 


QRASlllO 

c 

CHECK FOR NO WEIGHTING 

QRAS1120 

c 


GRAS1130 


10 IFCWTCID .LT. O.OD GO TO 80 

QRAS1140 

c 


QRAS1150 

c 

CHECK FOR ILLEGAL WEIGHTS 

QRAS1160 

c 


ORAS1170 


DO 20 I - 2, M 

QRAS1180 


IFCWTCID .LE. O.OD GO TO 30 

aRAS1190 


20 CONTINUE 

QRAS1200 


GO TO 40 

QRAS1210 


30 I ERR ' 3 

GRAS1220 


GO TO 160 

QRAS1230 

c 


QRAS1240 

c 

WEIGHT THE A AND B ARRAYS BY THE SQUARE ROOT 

GRAS1250 



C OF THE WEIGHT ARRAY. 

C 

40 DO 70 I - 1 , M 

WTCn « CQRTCWTCnD 
DO 50 J - 1, N 

ACI , J) - WTCU«ACI , J3 
50 CONTINUE 

DO 60 J - 1, IP 

ecl , j:) - wrcn^KBCi , jd 

60 CONTINUE 

70 CONTINUE 
80 CONTINUE 
C 

C CALL SQR0C2 TO DECOMPOSE A 

C 

CALL S0RDC2CA,MAXM,M,N,WKD 
C 

C CALL SGPEL2 TO SOLVE FOR IP RIGHT HAND SIDES 

C 

CALL SQRSL2<A,MAXM,M,N,MAXN, 1P,HK,E,X,RSD, JOB, 1ERP> 
IFCIERR .EQ. 0) GO TO 90 
lERR -= 2 
GO TO ISO 
90 CONTINUE 
C 

C COMPUTE THE SUM OF WEIGHTED SQUARES OF RESIDUALS. 

C 

1FCJ03 .EQ. n GO TO 140 
DO 110 J - 1 , IP 
SUMCJ5 - O.C 
DO iOO I - 1 , M 

smcjj ♦ R3DCL,JH:RSD(1>J> 


QRAS126G 

ORA51270 18 

QRAS1280 

QRAS1290 

QRAS1300 

QRAS1310 

QRAS1320 

QRAS1330 

QRAS1340 

QRAS1350 

QRAS1360 

QRAS1370 

GRA51380 

QRAS1390 

QRAS1400 

QRAS1410 

QRA51420 

QRAS1430 

QRAS1440 

QRAS1450 

QRAS14S0 

GRAS1470 

QRA51480 

QRAS1490 

QRAS1500 

QRAS1510 

GRA51520 

QRAS1530 

QRA51540 

QRAS1550 

QRAS1560 

QRAS157P 



iOO 


CONTI HUE 


QRAS1580 


110 CONTINUE 
C 

C COMPUTE UNWEIGHTED RESIDUALS 

C 

IFCWTCn .LT. 0.05 GO TO 160 
DO 130 1 - 1, M 

DO 120 J - 1, IP 

RSDC1,J5 - RSDCl,J5/WTCn 
120 CONTINUE 

130 CONTINUE 
140 .CONTINUE 

irCWTClJ .LT. 0.05 GO TO 160 
DO 150 I-1,M 

MTCI5 - WTcnnsTcn 
150 CONTINUE 
160 CONTINUE 
RETURN 
END 

SUBROUTINE S0.PPC2CX , LDX ,H ,P , QRAUX5 
INTEGER LDX,H,P 
REAL X(LDX,15 ,QRAUXa5 
C 

C SQRDC2 USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QP 

C FACTORIZATION OF AH N BY P MATRIX X. 

C 

C OH ENTRY 

C 

c X REALCLDX,P5. WHERE LDX .GE. N. 

C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE 

C COMPUTED. 


QRAS1590 

QRAS1600 

QRA51610 

QRA51620 

QRAS1630 

QRAS1640 

QRAS1650 

QRAS1660 

QRAS1670 

QRA51680 

QRA&1690 

QRAS1700 

QRAS1710 

QRAS1720 

QRAS1730 

QRAS1740 

QRAS1750 

QRAS17G0 

QRAS1770 

QRA3178C 

QRAS1790 

QRAS18CO 

QRAS1810 

CRAS1820 

QRAS1830 

QRAS1840 

QRAS1850 

QRAS1860 

QRAS1870 

QRAS18B0 



c 


QRA51890 


LDX INTEGER. 


QRAS1900 20 


LDX IS THE LEADING DIMENSION or THE ARRAY X, 


GRAS1310 


QRAS1920 


INTEGER. 


GRAS1930 


N IS THE HUMBER OF ROWS OF THE MATRIX X. 


QRAS1940 


QRAS1950 


INTEGER. 


QRAS19GO 


P IS THE NUMBER OF COLUMNS OF THE MATRIX X. 


QRAS1970 


QRA31980 


QRAS1990 


ON RETURN 


QRAS2000 


QRAS20I0 


X CONTAINS IN ITS UPPER TRIANGLE THE UPPER 


QRAS2020 


TRIANGULAR MATRIX R OF THE QR FACTOR IZAT1 OH. 


QRA52030 


BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM 


QRAS2040 


WHICH THE ORTHOGONAL PART OF THE DECOMPOSITION 


QRAS2050 


CAN BE RECOVERED. 


QRAS20E0 


GRAS2070 


QRAUV PEALCP^. 


QRAS2080 


QRAUX CONTAINS FURTHER ^’FORMATION REQUIRED TO RECOVER QRAS2090 


THE ORTHOGONAL PART OF THE DECOMPOSITION. 


QRAS2100 


QRAS2110 


0RA32120 


LINPACK SUBROUTINE SQRDC VERSION DATED 07/14/77, REVISED BY 


GRAS2130 


COMPUTER SCIENCES CORPORATION, HAMPTON, VA . 10/10/78. 


QRAS2140 


GRA52150 


BLAS SAXPYl ,SD0T1,SSCAL LRC NORMS 


GRA32160 


FORTRAN ABSSIGN, SORT, MOD 


QRAS2170 


QRAS2180 


INTERNAL VARIABLES 


ORAS2130 


0 




QRAS2200 




1MTE<5ER J,L,LP1 

QRAS2210 



21 


REAL SDOT.HRMXL.T 

QRAS2220 

c 


QRAS223D 

c 


QRAS22<40 

c 


QRAS2250 

c 

PERFORM the HOUEEHOLDEP REDUCTION OF X. 

QRAS2260 

c 


QRAS2270 


DO 190 L - 1, P 

QRAS2280 


QRAUXCD - O.CEO 

QRAS2290 


IF CL .EG. H> GO TO 170 

QRAS2300 

c 


QRAS2310 

c 

COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. 

QRAS2320 

c 


QRAS2330 


HLEH • N-L+1 

QRAS2340 


CALL NORMSCNLEM ,HLEN , 1 ,XCL ,L7 , 2 ,HRMXL3 

QRAS2350 


IF CHRMXL .EQ. O.OEO> C-0 TO 160 

QRAS2360 


IF CXCL.LD .HE. 0.0E9> NRMXL - SIGN CHRMXL , X CL ,L) 7 

QRAS2370 


CALL SSCALCN-L*l,1.0E0/NRMXL,XfL,L5 ,15 

GRAS2380 


XCL.L5 - l.OEO ♦ XCL,LD 

GRAS2390 

c 


ORAS24CO 

c 

APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS. 

QRAS2410 

c 


QRAS2420 


LPl - L ♦ 1 

QRAS2430 


IF CP .LT. LP15 GO TO 150 

QRAS2440 


DO 140 J - LPl, P 

QRAS2450 


T - -SDOTlCH-L+1 ,XCL,L5 ,XCL,J55/XCL,L5 

QRAS2460 


CALL SAXPYlCN-L + 1 ,T,XCL,L5 ,XCL, J5 5 

QRAS2470 


140 CONTINUE 

QR.AE2480 


150 CONTINUE 

QRAS2490 

c 


QRAS2500 


c 


QRAS2510 



C SAVE THE TRAHSrORMATION. 

GRAUXCL) - Xf.L,LD 
XCL,L> - -MRMXL 
160 CONTINUE 

170 CONTINUE 
180 CONTINUE 
RETURN 
END 


SUBROUTINE SQRSL2CK ,LDX ,N ,K ,LDB, IP ,QRAUX,Y ,BETA ,RSD , JOB, INFO^ 


INTEGER LDX,N,K,LDB,IP,JOD,INFO 

REAL XCLOX,n ,QRAUXCi: , YCLDX , n , BETACLD0 , 13 ,RSDCLDX,1D 


C 

c 

C 

c 

C 

C 

c 

C 

C 

c 

C 

C 

C 

c 

C 

C 

C 

C 

c 

c 

0 


SQRSL2 APPLIES THE OUTPUT OF THE SUBROUTINE SQRDC2 TO 
COMPUTE A SET OF IP LEAST SOUAPES SOLUTIONS AND RESIDUAL 
OUTPUT OF SQRDC2 IS THE DECOMPOSITION OF THE N BY X MATR 
X IN THE FORM 

X « Q M CR3 
(03 

MHERE Q IS ORTHOGONAL AND R IS UPPER TRIANGULAR. THIS 
INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAY X 
AMD THE ARRAY QRAUX. 

ON ENTRY 

X REALCLDX,F3, liHERE LDX .GE. N. 

X CONTAINS THE OUTPUT FROM SQRDC. 

LDX INTEGER. 

LDX IS THE LEADING DIMENSION OF THE ARRAY X. 


IX 


THE 


QRAS2520 

QRA52530 22 

QRAS2540 

QRAS2550 

QRAS2560 

QRAS2570 

QRAS2580 

QRAS2590 

QRAS2600 

QRAS2810 

QRAS2620 

QRAS2630 

QRAS2640 

ORAS2650 

QRAS2660 

QRAS2670 

QRAS2680 

QRAS269C 

QRAS2700 

CRAS2710 

QRAS2720 

QRAS2730 

QRA5274C 

QRAS2750 

QRAS2760 

QRAS2770 

QRAS2780 

GRAS2730 

QRAS2800 

QRAS2810 

QRAS2S20 


C.RAS2830 



c 

•H 

INTEGER. 

QRAS2840 

c 


H IS THE HUMBER OF ROWS OF .HE MATRIX X. 

QPA52850 

c 



aPf^S29S0 

c 

K 

INTEGER. 

QRA52870 

c 


K IS THE HUMBER OF COLUMNS OF THE MATRIX X. 

QRAS2880 

c 



QRA52890 

c 

LDB 

INTEGER. 

QRAS2900 

c 


LDB IS THE LEADING DIMENSION OF THE ARRAY BETA. 

QRAS2310 

c 



QRAS2920 

c 

IP 

INTEGER. 

QRAS2930 

c 


IP IS THE NUMBER OF RIGHT HAND SIDES. 

QRAS2940 

c 



QRAS2350 

c 

QRAUX 

REALCK) 

QRAS23S0 

0 


QRAUX CONTAINS THE OUTPUT FROM SQRDC2. 

QRAS2970 

c 



QRAS2980 

c 

Y 

REALCLDX.IP). 

QRAS2390 

c 


Y IS THE N BY IP RIGHT HAND SIDE MATRIX THAT IS 

QRA93000 

c 


MANIPULATED BY SQREL2. 

GRAS3010 

c 



QRAS3020 

c 



QRAS3030 

c 

JOB 

1 HTEGER . 

QRAS3040 

c 


JOB IS A PARAMETER THAT CONTROLS WHAT IS TO BE 

QRAS3050 

c 


COMPUTED. 

GRAS3060 

c 



QRAS3070 

c 


IF JOB .EG. 1 COMPUTE SOLUTICHS ONLY. 

QRAS3080 

c 


IF JOB .EQ. 2 COMPUTE RESIDUALS OHLY. 

QRAS3030 

c 


IF JOB .EQ. 3 COMPUTE SOLUTIONS AND RESIDUALS. 

QRAS3100 

c 



GRASSllO 

c 

ON RETURN 


QRAS3120 

c 



QRAS3130 

c 

BETA 

REALCLD0, IP). 

QRA33140 



c 

SETA COMTAIHS THE SOLUTIONS OF THE LEAST SQUARES 


QRAS3150 

c 

PROBLEMS 



QRAS3160 24 

c 

MINIMIZE HORMECYCn - XKBETACn:-, 1-1,2,. 

. • , IP 

QRAS3170 

c 

IF THEIR COMPUTATION HAS OEEN REQUESTED. 



QRAS3180 

c 




QRAS3190 

c 

RSO REALCLDX,IP5 



QRA93200 

c 

RSO CONTAINS THE LEAST SQUARES RESIDUALS 



QRAS3210 

c 

YCn - X*BETACn, 1-1,2,. ...IP 



QRAS5220 

c 

IF THEIR COMPUTATION HAS BEEN REQUESTED. 



QRAS3230 

c 




QRAS3240 

c 

INFO INTEGER 



QRAS3250 

c 

INFO IS ZERO UNLESS THE CALCULATION OF BETA 

HAS 

BEEM 

QRA33260 

c 

REQUESTED AND P IS SINGULAR, IN '4HICH CASE INFO 

IS 

QRAS3270 

c 

THE INDEX OF THE ^IRST ZERO DIAGONAL ELEMENT 

GF 

R. 

QRAS3290 

c 

IN THIS CASE BETA IS UNALTERED. 



QRA33290 

c 




QRA53300 

c 

LIHPACK SUBROUTII'E SQRSL VERSION DATED 07/14/77, REVISED 

DY 


QRA53310 

c 

COMPUTER SCIENCES CORPORATION, HAMPTON. VA. 10/10/73. 



QRAS3320 

c 




QRA53330 

c 

BLAS SAXPYl ,SCOPY,SDCTi 



QRAS3340 

c 

FORTRAN ABS.MINO.MOD 



QRAS3350 

c 




GR A333B0 

c 

INTERNAL VARIABLES 



QRAS3370 

c 




QRAS3380 


INTEGER I ,J, JJ.JU.KPl 



ORAS3390 


REAL SDOT.T.TEMP 



QRA53400 

c 




QRAS341D 

c 




QRAS3420 

c 

SET INFO FLAG 



GP>AS3430 

c 




ORAS3440 


INFO - 0 



GRAS3450 


jy_- MlNO(K,M-n 



QRAG3490 



QRAS3470 



SPECIAL ACTION WHEN H-l 

GRAS3480 



QRAS3490 


IF CJU .HE. 0) GO TO 20 

QRAS3500 


IFCXCl,!^ .HE. O.C) GO TO 5 

QRAS3510 


INFO - 1 

QRA53520 


GO TO 220 

QRAS3530 

5 

CONTINUE 

QRAS3540 


DO 10 L • 1, IP 

QRAS3550 


IFCJ03 .NE. 2) BETACl.L) - YC1,L1/XC1,1) 

GRAS3560 


IFCJOB .NE. 10 RS0C1,L5 - O.OEO 

QRA33570 

10 

CONTINUE 

QRAS3580 


GO TO 220 

QRAS3590 

20 

CONTINUE 

QRA33600 



QRA53610 


COMPUTE TRANS CQ5*Y 

QRAS3620 



QRAS363G 


DO 50 J - 1, JU 

0.RAS3640 


IF CQRAUXCJ5 .EQ. 0.CE05 GO TO 40 

QRASoSSO 


TEMP - XCJ.J) 

GRA53660 


XCJ.J) - QRAUXCJ5 

QRA53S70 


DO 30 L - 1, IP 

QRAS3680 


T - -SD0T1CH-J«1,XCJ,J5 ,YCJ,L))/XCJ,JD 

GRAS36O0 


CALL SAXPVlCH-J+1 ,T,XCJ,J) ,YCJ,0) 

QRAS3700 

30 

CONTINUE 

QR AS3 / 1 C 


XCJ.Jl - TEMP 

GRAS3720 

40 

CONTINUE 

QRAS3730 

50 

CONTINUE 

QRAS3740 


KPl - K + 1 

QRAS3750 


IF (JOB .EO. 1 .OR. K -EQ. H) GO 70 70 

QRAS37GO 


DO GO L 


1, IP 


GPAS3770 



CALI iCOFYCH-K.YCKPl ,Li , 1 .RSDCKPl ,L) ,n 


QRAS3780 


60 

COHTIHUE 


QRAS3790 

70 

CONTINUE 


QRAS3800 


IF (JOB .EQ. 2!) GO TO 120 


QRAS3P1D 

C 



QRAS3820 

C 

COMPUTE BETA 


QRAS3830 

C 



QRA93B40 


DO 75 L - 1, IP 


QRA53850 


CALL SC0PY(K,YC1 ,L> ,1,BETAC1 >L) ,1} 


QRAS3860 

75 

CONTINUE 


QRAS3870 


DO 100 JJ - 1, K 


QRAS3880 


J • K - JJ * 1 


QRAS3890 


IF CXCJ.J) .ME. O.OEO^ GO TO 80 


aRAS3900 


INFO - J 


QRA33910 

C 

EXIT 


QRAS3920 


GO TO 220 


QRAS3930 

80 

CONTINUE 


QRAS3940 


DO 95 L - 1, IP 


QRAS3950 


BETACJ,L) - BETACJ,LD /XCJ, J3 


QRAS3960 


IF CJ -EQ. 1) GO TO 90 


QRAS3970 


T « -BETACJ>L3 


GRAS3980 


CALL SAXPYIC J-1>T,XC1,J) ,BETAC1 ,L0D 


QPAS39S0 

90 

COHTIHUE 


GRAS4000 

95 

COHTIHUE 


QRA54D10 

100 

CONTINUE 


GRAS^020 

110 

CONTINUE 


QRAS4030 

120 

COHTIHUE 


QRAS4040 


IF (JOB .EG. n GO TO 210 


QRAS4050 

c 



QRAS40G0 

c 

COMPUTE PSD IF REQUIRED 


QRAS4070 

c 



QRAS4080 


DO IGO L - 1 , IP 


QRAS4090 



150 

IGO 


170 

130 

200 

210 

220 


Du 150 I - i, K 
RSDCIjLj ■ o.oeo 
COHTIKuE 
COHTIHiiE 

Du 200 JJ - 1, JU 
J - JU - JJ ^ i 

IF CQRmUaCJ^ .EQ. O.GEOO GO TO ISO 
TEMP * XCJ,Jj 
aCJ,Jj - GRhUXCJj 
DO 170 L - 1, IP 

T » "SDOTl CH- J-** 1 ^X C J , JO , RSD C J ,L5 j /X C J , Jj 
CALL SHXFYiCri- J^i ,T,XCJ, JO ,RSDC J,LOO 
COKT I HUE 
XCviyJO • TEMP 
CONTI HUE 
COHT 1 HUE 
COKT I HUE 
COKT i HUE 
RETURH 
END 


QRASmOG 

GRASmiO 

QRh34120 

QRASmoO 

QRAS-imG 

QRAS4i50 

QRAS416G 

QRA34170 

QRA54 180 

GRA3413G 

GRA54200 

QRA5421G 

QRA5422G 

QRA3423G 

uRA54240 

QRAS425D 

QRA542G0 

QRA34270 

QRA34260 

URA34290 
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